##### Load Data

	cd(datadir)

	# Packages
	using DataFrames
	using StatsBase
	using Optim
	using Calculus
	using Distributions
	using Random
	using DelimitedFiles
	using CSV
	using LinearAlgebra
	
	# Load Data 
	if network_flag
		if fill_missing
			data = CSV.read("ca_julia_data_JUN182021_net_fillmissing.csv",DataFrame, header=true); # this has 4 fields (choice, premium, plan name, and household id)
			households = CSV.read("ca_household_characteristics_JUN182021_net_fillmissing.csv",DataFrame, header = true) ; # household characteristics 
			plans = CSV.read("ca_plan_characteristics_JUN182021_net.csv",DataFrame, header = true) ; # plan characteristics (does not include same plans in different rating areas and years)
		elseif keep_missing
			data = CSV.read("ca_julia_data_JUN182021_net_keepmissing.csv",DataFrame, header=true); # this has 4 fields (choice, premium, plan name, and household id)
			households = CSV.read("ca_household_characteristics_JUN182021_net_keepmissing.csv",DataFrame, header = true) ; # household characteristics 
			plans = CSV.read("ca_plan_characteristics_JUN182021_net.csv",DataFrame, header = true) ; # plan characteristics (does not include same plans in different rating areas and years)
		else
			data = CSV.read("ca_julia_data_JUN182021_net.csv",DataFrame, header=true); # this has 4 fields (choice, premium, plan name, and household id)
			households = CSV.read("ca_household_characteristics_JUN182021_net.csv",DataFrame, header = true) ; # household characteristics 
			plans = CSV.read("ca_plan_characteristics_JUN182021_net.csv",DataFrame, header = true) ; # plan characteristics (does not include same plans in different rating areas and years)
		end
	else
		data = CSV.read("ca_julia_data_AUG032019_small.csv",DataFrame,header=true); # this has 4 fields (choice, premium, plan name, and household id)
		households = CSV.read("ca_household_characteristics_AUG032019_small.csv",DataFrame,header = true) ; # household characteristics 
		plans = CSV.read("ca_plan_characteristics_AUG032019_small.csv",DataFrame,header = true) ; # plan characteristics (does not include same plans in different rating areas and years)
	end	
	Random.seed!(6);

	# Output file names
	
	if nested
		filename = "nested"; # file to save results
		choice_output_filename = "nested_choice_probs";
		elasticities_filename = "nested_elasticities";
		inertia_output_filename = "nested_switching costs";
		avg_param_output_filname = "nested_average_param_output";
	else
		filename = "logit"; # file to save results
		choice_output_filename = "logit_choice_probs";
		elasticities_filename = "logit_elasticities";
		inertia_output_filename = "logit_switching costs";
		avg_param_output_filname = "logit_average_param_output";
	end
	
	if eq_robust_flag
		filename = string(filename,"_eqrobust");
		choice_output_filename = string(choice_output_filename,"_eqrobust");
		elasticities_filename = string(elasticities_filename,"_eqrobust");
		inertia_output_filename = string(inertia_output_filename,"_eqrobust");
		avg_param_output_filname = string(avg_param_output_filname,"_eqrobust");
	end
	
	if inattention
		filename = string(filename,"_inattention");
		choice_output_filename = string(choice_output_filename,"_inattention");
		elasticities_filename = string(elasticities_filename,"_inattention");
		inertia_output_filename = string(inertia_output_filename,"_inattention");
		avg_param_output_filname = string(avg_param_output_filname,"_inattention");
	end
	
	if sequence_flag
		filename = string(filename,"_seq");
		choice_output_filename = string(choice_output_filename,"_seq");
		elasticities_filename = string(elasticities_filename,"_seq");
		inertia_output_filename = string(inertia_output_filename,"_seq");
		avg_param_output_filname = string(avg_param_output_filname,"_seq");
	end
	
	if random_parameters
		filename = string(filename,"_random");
		choice_output_filename = string(choice_output_filename,"_random");
		elasticities_filename = string(elasticities_filename,"_random");
		inertia_output_filename = string(inertia_output_filename,"_random");
		avg_param_output_filname = string(avg_param_output_filname,"_random");
	end
	
	if control_function
		filename = string(filename,"_cf");
		choice_output_filename = string(choice_output_filename,"_cf");
		elasticities_filename = string(elasticities_filename,"_cf");
		inertia_output_filename = string(inertia_output_filename,"_cf");
		avg_param_output_filname = string(avg_param_output_filname,"_cf");
	end
	
	if add_av_interactions
		filename = string(filename,"_av");
		choice_output_filename = string(choice_output_filename,"_av");
		elasticities_filename = string(elasticities_filename,"_av");
		inertia_output_filename = string(inertia_output_filename,"_av");
		avg_param_output_filname = string(avg_param_output_filname,"_av");
	end	
	
	if add_insurer_market_fe
		filename = string(filename,"_imfe");
		choice_output_filename = string(choice_output_filename,"_imfe");
		elasticities_filename = string(elasticities_filename,"_imfe");
		inertia_output_filename = string(inertia_output_filename,"_imfe");
		avg_param_output_filname = string(avg_param_output_filname,"_imfe");
	end	
	
	if add_silver
		filename = string(filename,"_silver");
		choice_output_filename = string(choice_output_filename,"_silver");
		elasticities_filename = string(elasticities_filename,"_silver");
		inertia_output_filename = string(inertia_output_filename,"_silver");
		avg_param_output_filname = string(avg_param_output_filname,"_silver");
	end	
	
	if network_flag
		filename = string(filename,"_netpair");
		choice_output_filename = string(choice_output_filename,"_netpair");
		elasticities_filename = string(elasticities_filename,"_netpair");
		inertia_output_filename = string(inertia_output_filename,"_netpair");
		avg_param_output_filname = string(avg_param_output_filname,"_netpair");
	end
	
	if keep_years_in_panel > 0
		filename = string(filename,"_",keep_years_in_panel);
		choice_output_filename = string(choice_output_filename,"_",keep_years_in_panel);
		elasticities_filename = string(elasticities_filename,"_",keep_years_in_panel);
		inertia_output_filename = string(inertia_output_filename,"_",keep_years_in_panel);
		avg_param_output_filname = string(avg_param_output_filname,"_",keep_years_in_panel);
	end	
	
	if use_newton
		filename = string(filename,"_newton");
		choice_output_filename = string(choice_output_filename,"_newton");
		elasticities_filename = string(elasticities_filename,"_newton");
		inertia_output_filename = string(inertia_output_filename,"_newton");
		avg_param_output_filname = string(avg_param_output_filname,"_newton");
	end
	
	if smooth_flag
		filename = string(filename,tau);
	end
	
	filename = string(filename,".csv");
	choice_output_filename = string(choice_output_filename,".csv");
	elasticities_filename = string(elasticities_filename,".csv");
	inertia_output_filename = string(inertia_output_filename,".csv");
	avg_param_output_filname = string(avg_param_output_filname,".csv");
		
	# Reduce Sample
	if eq_robust_flag
		keep_households = intersect(findall(convert(Array,households[!,:year]) .>= 2016),
			findall(convert(Array,households[!,:year]) .<= year_to_run));
	else
		keep_households = findall(convert(Array,households[!,:year]) .<= year_to_run);
	end
	if rich_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:FPL]) .>= 2));
	end
	if poor_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:FPL]) .< 2));
	end
	if old_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:mean_age]) .>= 45));
	end
	if young_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:mean_age]) .< 45));
	end
	if lapse_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:ever_lapsed]) .== 1));
	end
	if no_lapse_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:ever_lapsed]) .== 0));
	end
	if churn_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:ever_churned]) .== 1));
	end
	if no_churn_flag
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:ever_churned]) .== 0));
	end
	
	if keep_years_in_panel > 0
		keep_households = intersect(keep_households,findall(convert(Array,households[!,:years_in_panel_1418]) .== keep_years_in_panel));
	end
	
	if keep_years_in_panel == 1
		add_previous_demographis_interactions = false;
		add_previous_firm_interactions = false;
	end
	
	# Take smaller sample
	if shrink_sample
		num_draws = Int64(round(sample_size*length(unique(households[keep_households,:household_id]))));
		keep_sequence_ids = sort(sample(unique(households[keep_households,:household_id]),num_draws,replace=false));
		keep_households = intersect(keep_households,findall(in(keep_sequence_ids),households[!,:household_id]));
	end
	
	if keep_two_choices | keep_three_choices
		years_in_panel = combine(groupby(households[keep_households,:],:household_id), :household_id => length)[!,:household_id_length];
		time_in_panel = zeros(size(households)[1]);
		unique_hseq_ids = unique(households[keep_households,:household_id]);
		for s in 1:length(unique_hseq_ids)
			time_in_panel[findall(households[!,:household_id] .== unique_hseq_ids[s])] .= years_in_panel[s];
		end
		
		if keep_two_choices
			# Delete households only in panel for 1 year
			keep_households = intersect(keep_households,findall(time_in_panel .> 1));
			
			# Delete households' choices in years 3 and above
			first_year_in_panel = combine(groupby(households[keep_households,:],:household_id), :year => minimum)[!,:year_minimum];
			in_first_two_years = zeros(size(households)[1]);
			unique_hseq_ids = unique(households[keep_households,:household_id]);
			for s in 1:length(unique_hseq_ids)
				in_first_two_years[intersect(findall(households[!,:household_id] .== unique_hseq_ids[s]),findall(households[!,:year] .<= first_year_in_panel[s] + 1))] .= 1;
			end
			keep_households = intersect(keep_households,findall(in_first_two_years .== 1));
			
			# Some Households are in sample for two non-consecutive years - delete them
			years_in_panel = combine(groupby(households[keep_households,:],:household_id), :household_id => length)[!,:household_id_length];
			time_in_panel = zeros(size(households)[1]);
			unique_hseq_ids = unique(households[keep_households,:household_id]);
			for s in 1:length(unique_hseq_ids)
				time_in_panel[findall(households[!,:household_id] .== unique_hseq_ids[s])] .= years_in_panel[s];
			end
			keep_households = intersect(keep_households,findall(time_in_panel .> 1));
		end
		
		if keep_three_choices
			# Delete households only in panel for 1 or 2 years
			keep_households = intersect(keep_households,findall(time_in_panel .> 2));
			
			# Delete households' choices in years 4 and above
			first_year_in_panel = combine(groupby(households[keep_households,:],:household_id), :year => minimum)[!,:year_minimum];
			in_first_three_years = zeros(size(households)[1]);
			unique_hseq_ids = unique(households[keep_households,:household_id]);
			for s in 1:length(unique_hseq_ids)
				in_first_three_years[intersect(findall(households[!,:household_id] .== unique_hseq_ids[s]),findall(households[!,:year] .<= first_year_in_panel[s] + 2))] .= 1;
			end
			keep_households = intersect(keep_households,findall(in_first_three_years .== 1));
			
			# Some Households are in sample for three non-consecutive years - delete them
			years_in_panel = combine(groupby(households[keep_households,:],:household_id), :household_id => length)[!,:household_id_length];
			time_in_panel = zeros(size(households)[1]);
			unique_hseq_ids = unique(households[keep_households,:household_id]);
			for s in 1:length(unique_hseq_ids)
				time_in_panel[findall(households[!,:household_id] .== unique_hseq_ids[s])] .= years_in_panel[s];
			end
			keep_households = intersect(keep_households,findall(time_in_panel .> 2));
		end
	end
	
	# Reassign household numbers
	new_household_numbers = zeros(Int64,size(households)[1]);
	new_household_numbers[keep_households] = collect(1:1:length(keep_households));
	data[!,:household_number] = new_household_numbers[convert(Array,data[!,:household_number])];
	
	# Delete records outside sample
	writedlm("data_1percsample.csv",findall(data[!,:household_number] .> 0))
	data = data[data[!,:household_number] .> 0,:];
	households = households[keep_households,:];
	
	
	subsidies = convert(Array,data[!,:subsidy]);
	penalties = convert(Array,data[!,:penalty]);
	
	# Specify Covariates (NOTE: control function approach should include "residual","eta"; random parameters should include random parameter covariates)
	all_covariates = convert(Array,starting_point[!,:Variable];);
	if nested
		pop!(all_covariates);
	end
	if inattention
		premium_shock_index = findall(all_covariates .== "premium_shock")[1];
		search_covariates = all_covariates[premium_shock_index:length(all_covariates)];
		all_covariates = all_covariates[1:(premium_shock_index-1)];
	end
	
	if add_complex_interactions & add_silver
		product_covariates = ["Premium","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];	
	elseif add_complex_interactions
		product_covariates = ["Premium","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];	
	elseif add_silver
		product_covariates = ["Premium","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net","Anthem_HMO"];	
	else
		product_covariates = ["Premium","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net","Anthem_HMO"];	
	end
	
	if random_parameters
		if inattention & add_silver
			random_variables = ["Anthem","Blue_Shield","Kaiser","Health_Net","AV","Silver"];
			random_parameter_product_indices = findall(in(random_variables),product_covariates);
		elseif inattention & add_av_interactions
			random_variables = ["Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = findall(in(random_variables),product_covariates);
		elseif inattention
			random_variables = ["Anthem","Blue_Shield","Kaiser","Health_Net","AV"];
			random_parameter_product_indices = findall(in(random_variables),product_covariates);
		elseif add_complex_interactions & add_av_interactions
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		elseif add_complex_interactions & add_silver
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net","AV","Silver"];
			#random_covariance_variables = ["Anthem_random","BS_Anthem_cov","Kaiser_Anthem_cov","HN_Anthem_cov","AV_Anthem_cov",
			#								"Blue_Shield_random","Kaiser_BS_cov","HN_BS_cov","AV_BS_cov","Kaiser_random","HN_Kaiser_cov",
			#								"AV_Kaiser_cov","Health_Net_random","AV_HN_cov","AV_random"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		elseif add_complex_interactions 
			random_variables = ["intercept","Anthem","Blue_Shield","Kaiser","Health_Net","AV"];
			random_covariance_variables = ["Anthem_random","BS_Anthem_cov","Kaiser_Anthem_cov","HN_Anthem_cov","AV_Anthem_cov",
											"Blue_Shield_random","Kaiser_BS_cov","HN_BS_cov","AV_BS_cov","Kaiser_random","HN_Kaiser_cov",
											"AV_Kaiser_cov","Health_Net_random","AV_HN_cov","AV_random"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		elseif add_silver
			random_variables = ["intercept","AV","Silver","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		else
			random_variables = ["intercept","AV","HMO","Anthem","Blue_Shield","Kaiser","Health_Net"];
			random_parameter_product_indices = cat(dims=1,findall(all_covariates .== "intercept")[1],findall(in(random_variables),product_covariates));
		end
		random_product_covariates = string.(random_variables,"_random");
		
		random_parameter_indices = findall(in(random_product_covariates),all_covariates);
		if length(random_parameter_indices) == 1
			random_parameter_indices = random_parameter_indices[1];
		end
	end
	
	# Estimation variables and parameters
		
		choices = convert(Array,data[!,:choice]);
		data_weights = convert(Array,data[!,:weight]);
		uninsured_plans = convert(Array,data[!,:uninsured_plan]);
		household_numbers = convert(Array,data[!,:household_number]);	
		
		if inattention
			entrant_weights = data_weights .* households[household_numbers,:entrant];
			switcher_weights = data_weights .* households[household_numbers,:switcher];
			nonswitcher_weights = data_weights .* households[household_numbers,:nonswitcher];
		end
		
		unique_plans = unique(data[!,:plan_name]);	
		plan_indicator = trues(size(plans)[1]);
		for i in eachindex(plan_indicator)
			plan_indicator[i] = in(plans[i,:Plan_Name],unique_plans);
		end
		plans = plans[findall(plan_indicator),:];
			
		# Indicator for nested logit gradient and hessian
		I_ex_nest_household = 1 .- choices[uninsured_plans .== 1];
		I_ex_nest = I_ex_nest_household[household_numbers];	
			
		# Dataframe for determining choice sequence probabilities
		observed_choice_indices = findall(choices .== 1);
		observed_choice_household_ids = households[household_numbers[observed_choice_indices],:household_id];	
		hseq_probs = DataFrame();
		hseq_probs[!,:household_id] = observed_choice_household_ids;
		hseq_probs[!,:prob] = zeros(length(observed_choice_household_ids));
		hseq_probs[!,:weight] = households[household_numbers[observed_choice_indices],:weight];	
		hseq_probs[!,:X] = zeros(length(observed_choice_household_ids));
		#hseq_weights = combine(groupby(hseq_probs,:household_id), :weight => minimum)[!,:weight_minimum];
		#hseq_weights = combine(groupby(hseq_probs,:household_id), :weight => sum)[!,:weight_sum];
		hseq_weights = combine(groupby(hseq_probs,:household_id), :weight => mean)[!,:weight_mean];
		
		if sequence_exp_flag 
			household_hseq_weights = zeros(size(households)[1]);
			unique_hseq_ids = unique(hseq_probs[!,:household_id]);
			for s in 1:length(unique_hseq_ids)
				household_hseq_weights[findall(hseq_probs[!,:household_id] .== unique_hseq_ids[s])] .= hseq_weights[s];
			end
			data_weights = household_hseq_weights[household_numbers];
		end
		
		
	# Form X covariate matrix
	function form_Xsearch_matrix(search_covariates)	

		Xsearch = zeros(Float64,size(households)[1],length(search_covariates));
		Xsearch[:,findall(search_covariates .== "premium_shock")[1]] = households[!,:premium_shock];
		Xsearch[:,findall(search_covariates .== "hh_250to400")[1]] = ((households[!,:FPL] .> 2.50) .& (households[!,:FPL] .<= 4)) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_gt400")[1]] = (households[!,:FPL] .> 4) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_0to17")[1]] = (households[!,:perc_0to17]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_18to34")[1]] = (households[!,:perc_18to25] .+ households[!,:perc_26to34]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_35to54")[1]] = (households[!,:perc_35to44] .+ households[!,:perc_45to54]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Black")[1]] = (households[!,:perc_black]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Hispanic")[1]] = (households[!,:perc_hispanic]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Asian")[1]] = (households[!,:perc_asian]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_Other")[1]] = (households[!,:perc_other]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_male")[1]] = (households[!,:perc_male]) .* 1; 
		Xsearch[:,findall(search_covariates .== "hh_family")[1]] = (households[!,:household_size] .> 1) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2015")[1]] = (households[!,:year] .== 2015) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2016")[1]] = (households[!,:year] .== 2016) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2017")[1]] = (households[!,:year] .== 2017) .* 1; 
		#Xsearch[:,findall(search_covariates .== "hh_2018")[1]] = (households[!,:year] .== 2018) .* 1; 
		Xsearch[:,findall(search_covariates .== "intercept_search")[1]] .= 1; 
		
		return(Xsearch)
	end
	
	# Form X covariate matrix
	function form_X_matrix(all_covariates)
	
		prem_ind = findall(all_covariates .== "Premium")[1];
		if !inattention & (keep_years_in_panel != 1)
			prev_ind = findall(all_covariates .== "previous_choice")[1];
			#silver_ind = findall(all_covariates .== "Silver")[1];
		end
		av_ind = findall(all_covariates .== "AV")[1];
		
		X = zeros(Float64,length(choices),length(all_covariates));
		X[:,prem_ind] = convert(Array,data[!,:premium]) ./ households[household_numbers,:enrollees];
		if !inattention & (keep_years_in_panel != 1)
			X[:,prev_ind] = convert(Array,data[!,:previous_choice]);
		end
		X[:,findall(all_covariates .== "AV")[1]] = convert(Array,data[!,:av]);
		X[:,findall(all_covariates .== "intercept")] = 1 .- uninsured_plans;		
		
		# Add product characteristics
		plans[!,:PPO] = plans[!,:HMO];
		for plan in plans[!,:Plan_Name]
			plan_indices = findall(data[!,:plan_name] .== plan);
			for covariate in setdiff(product_covariates,["Premium","AV","Anthem_HMO"])
				X[plan_indices,findall(all_covariates .== covariate)[1]] .= 
					plans[findall(plans[!,:Plan_Name] .== plan)[1],findall(names(plans) .== covariate)[1]];
			end
		end
		
		if !add_complex_interactions
			X[:,findall(all_covariates .== "Anthem_HMO")[1]] = X[:,findall(all_covariates .== "Anthem")[1]] .* X[:,findall(all_covariates .== "HMO")[1]];
		end
		
		# Add premium interactions
		if add_premium_interactions
			if (rich_flag | poor_flag | old_flag | young_flag)
				X[:,findall(all_covariates .== "Premium_male")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "Premium_family")[1]] = X[:,prem_ind] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "Premium_Asian")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "Premium_Black")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "Premium_Hispanic")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "Premium_Other")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_other]);
			else
				X[:,findall(all_covariates .== "Premium_250to400")[1]] = X[:,prem_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 
				X[:,findall(all_covariates .== "Premium_gt400")[1]] = X[:,prem_ind] .* (households[household_numbers,:FPL] .> 4) ; 
				if add_complex_interactions
					X[:,findall(all_covariates .== "Premium_0to34")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_0to17] .+ 
						households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				else
					X[:,findall(all_covariates .== "Premium_0to17")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_0to17]);
					X[:,findall(all_covariates .== "Premium_18to34")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				end
				X[:,findall(all_covariates .== "Premium_35to54")[1]] =  X[:,prem_ind] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "Premium_male")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "Premium_family")[1]] = X[:,prem_ind] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "Premium_Asian")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "Premium_Black")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "Premium_Hispanic")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "Premium_Other")[1]] = X[:,prem_ind] .* (households[household_numbers,:perc_other]);
			end
			
			if add_complex_interactions
				X[:,findall(all_covariates .== "Premium_250to400_0to34")[1]] = X[:,prem_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]); 
				X[:,findall(all_covariates .== "Premium_gt400_0to34")[1]] = X[:,prem_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);  
				X[:,findall(all_covariates .== "Premium_250to400_35to54")[1]] = X[:,prem_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]); 
				X[:,findall(all_covariates .== "Premium_gt400_35to54")[1]] = X[:,prem_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);  
			end
		end
	
		# AV interactions
		if add_av_interactions
			X[:,findall(all_covariates .== "AV_250to400")[1]] = X[:,av_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 
			X[:,findall(all_covariates .== "AV_gt400")[1]] = X[:,av_ind] .* (households[household_numbers,:FPL] .> 4) ; 
			X[:,findall(all_covariates .== "AV_0to17")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_0to17]);
			X[:,findall(all_covariates .== "AV_18to34")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
			X[:,findall(all_covariates .== "AV_35to54")[1]] =  X[:,av_ind] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
			X[:,findall(all_covariates .== "AV_male")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_male]);
			X[:,findall(all_covariates .== "AV_family")[1]] = X[:,av_ind] .* (households[household_numbers,:household_size] .> 1);
			X[:,findall(all_covariates .== "AV_Asian")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_asian]);
			X[:,findall(all_covariates .== "AV_Black")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_black]);
			X[:,findall(all_covariates .== "AV_Hispanic")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_hispanic]);
			X[:,findall(all_covariates .== "AV_Other")[1]] = X[:,av_ind] .* (households[household_numbers,:perc_other]);
			
			if add_complex_interactions
				X[:,findall(all_covariates .== "AV_250to400_0to34")[1]] = X[:,av_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]); 
				X[:,findall(all_covariates .== "AV_gt400_0to34")[1]] = X[:,av_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);  
				X[:,findall(all_covariates .== "AV_250to400_35to54")[1]] = X[:,av_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]); 
				X[:,findall(all_covariates .== "AV_gt400_35to54")[1]] = X[:,av_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);  
			end
		end
	
		# Add inertia interactions
		if add_previous_demographics_interactions & !inattention
			if (rich_flag | poor_flag | old_flag | young_flag)
				X[:,findall(all_covariates .== "previous_choice_male")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "previous_choice_family")[1]] = X[:,prev_ind] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "previous_choice_Asian")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "previous_choice_Black")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "previous_choice_Hispanic")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "previous_choice_Other")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_other]);
			else
				X[:,findall(all_covariates .== "previous_choice_250to400")[1]] = X[:,prev_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) ; 
				X[:,findall(all_covariates .== "previous_choice_gt400")[1]] = X[:,prev_ind] .* (households[household_numbers,:FPL] .> 4) ; 
				if add_complex_interactions
					X[:,findall(all_covariates .== "previous_choice_0to34")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_0to17] .+
						households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				else
					X[:,findall(all_covariates .== "previous_choice_0to17")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_0to17]);
					X[:,findall(all_covariates .== "previous_choice_18to34")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				end
				X[:,findall(all_covariates .== "previous_choice_35to54")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "previous_choice_male")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "previous_choice_family")[1]] = X[:,prev_ind] .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "previous_choice_Asian")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "previous_choice_Black")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_black]);
				X[:,findall(all_covariates .== "previous_choice_Hispanic")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "previous_choice_Other")[1]] = X[:,prev_ind] .* (households[household_numbers,:perc_other]);
			end
			
			if add_complex_interactions
				X[:,findall(all_covariates .== "previous_choice_250to400_0to34")[1]] = X[:,prev_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]); 
				X[:,findall(all_covariates .== "previous_choice_gt400_0to34")[1]] = X[:,prev_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_0to17] .+ households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);  
				X[:,findall(all_covariates .== "previous_choice_250to400_35to54")[1]] = X[:,prev_ind] .* ((households[household_numbers,:FPL] .> 2.50) .& (households[household_numbers,:FPL] .<= 4)) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]); 
				X[:,findall(all_covariates .== "previous_choice_gt400_35to54")[1]] = X[:,prev_ind] .* (households[household_numbers,:FPL] .> 4) .*
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);  
			end
		end
		
		if add_previous_firm_interactions & !inattention & !add_complex_interactions
			X[:,findall(all_covariates .== "previous_choice_Anthem")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "previous_choice_Blue_Shield")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "previous_choice_Kaiser")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "previous_choice_Health_Net")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Health_Net")[1]];			
			X[:,findall(all_covariates .== "previous_choice_HMO")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "HMO")[1]];
			X[:,findall(all_covariates .== "previous_choice_AV")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "AV")[1]];
			if add_silver
				X[:,findall(all_covariates .== "previous_choice_silver")[1]] = X[:,prev_ind] .* X[:,findall(all_covariates .== "Silver")[1]];
			end
		end
	
		if network_flag
			mean_breadth = sum(data[uninsured_plans .== 0,:network_breadth] .* data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight])/
				sum(data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight]);
			mean_inclus = sum(data[uninsured_plans .== 0,:network_inclus] .* data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight])/
				sum(data[uninsured_plans .== 0,:choice] .* data[uninsured_plans .== 0,:weight]);
				
			#X[:,findall(all_covariates .== "breadth")[1]] = convert(Array,data[!,:network_breadth]) .- 
			#	mean_breadth * (1 .- uninsured_plans);
			#X[:,findall(all_covariates .== "inclus")[1]] = convert(Array,data[!,:network_inclus]) .- 
			#	mean_inclus * (1 .- uninsured_plans);
			X[:,findall(all_covariates .== "pairinclus")[1]] = convert(Array,data[!,:network_pairinclus]);
			
			#X[:,findall(all_covariates .== "previous_choice_breadth")[1]] = 
			#	X[:,findall(all_covariates .== "previous_choice")[1]] .* 
			#		(convert(Array,data[!,:network_breadth]) .- mean_breadth * (1 .- uninsured_plans));
			#X[:,findall(all_covariates .== "previous_choice_inclus")[1]] = 
			#	X[:,findall(all_covariates .== "previous_choice")[1]] .* 
			#		(convert(Array,data[!,:network_inclus]) .- mean_inclus * (1 .- uninsured_plans));		
		end
	
		# Add intercepts
		if add_intercepts
			if add_complex_interactions
				if !eq_robust_flag
					X[:,findall(all_covariates .== "intercept_year2015")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2015);
					X[:,findall(all_covariates .== "intercept_year2016")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2016);
				end
				X[:,findall(all_covariates .== "intercept_year2017")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2017);
				X[:,findall(all_covariates .== "intercept_year2018")] = (1 .- uninsured_plans) .* (households[household_numbers,:year] .== 2018);
				#X[:,findall(all_covariates .== "intercept_0to34")] = (1 .- uninsured_plans)  .* (households[household_numbers,:perc_0to17] .+
				#	households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				#X[:,findall(all_covariates .== "intercept_35to54")] = (1 .- uninsured_plans) .* 
				#	(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "HMO_0to34")[1]] = X[:,findall(all_covariates .== "HMO")[1]] .* (households[household_numbers,:perc_0to17] .+
					households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);
				X[:,findall(all_covariates .== "HMO_35to54")[1]] = X[:,findall(all_covariates .== "HMO")[1]] .* 
					(households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
			else
				X[:,findall(all_covariates .== "intercept_male")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_male]);
				X[:,findall(all_covariates .== "intercept_family")] = (1 .- uninsured_plans) .* (households[household_numbers,:household_size] .> 1);
				X[:,findall(all_covariates .== "intercept_0to17")] = (1 .- uninsured_plans)  .* (households[household_numbers,:perc_0to17]);
				X[:,findall(all_covariates .== "intercept_18to34")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_18to25] .+ households[household_numbers,:perc_26to34]);	
				X[:,findall(all_covariates .== "intercept_35to54")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_35to44] .+ households[household_numbers,:perc_45to54]);
				X[:,findall(all_covariates .== "intercept_250to400")] = (1 .- uninsured_plans) .* ((households[household_numbers,:FPL] .> 2.5) .& (households[household_numbers,:FPL] .<= 4)) ;						
				X[:,findall(all_covariates .== "intercept_gt400")] = (1 .- uninsured_plans) .* (households[household_numbers,:FPL] .> 4) ; 
				X[:,findall(all_covariates .== "intercept_asian")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_asian]);
				X[:,findall(all_covariates .== "intercept_black")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_black]);	
				X[:,findall(all_covariates .== "intercept_hispanic")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_hispanic]);
				X[:,findall(all_covariates .== "intercept_other")] = (1 .- uninsured_plans) .* (households[household_numbers,:perc_other]);
			end
		end
		
		# Add rating areas
		if add_rating_areas
			X[:,findall(all_covariates .== "intercept_ra2")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2);
			X[:,findall(all_covariates .== "intercept_ra3")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3);
			X[:,findall(all_covariates .== "intercept_ra4")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4);
			X[:,findall(all_covariates .== "intercept_ra5")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5);
			X[:,findall(all_covariates .== "intercept_ra6")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 6);
			X[:,findall(all_covariates .== "intercept_ra7")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7);
			X[:,findall(all_covariates .== "intercept_ra8")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8);
			X[:,findall(all_covariates .== "intercept_ra9")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 9);
			X[:,findall(all_covariates .== "intercept_ra10")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10);
			X[:,findall(all_covariates .== "intercept_ra11")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 11);
			X[:,findall(all_covariates .== "intercept_ra12")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 12);
			X[:,findall(all_covariates .== "intercept_ra14")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14);
			X[:,findall(all_covariates .== "intercept_ra15")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15);
			X[:,findall(all_covariates .== "intercept_ra16")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16);
			X[:,findall(all_covariates .== "intercept_ra17")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17);
			X[:,findall(all_covariates .== "intercept_ra18")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18);
			X[:,findall(all_covariates .== "intercept_ra19")[1]] = (1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19);
		end	
		
		# Add insurer-market fixed effects
		if add_insurer_market_fe
			X[:,findall(all_covariates .== "intercept_ra2_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra2_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra2_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra2_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 2) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];

			X[:,findall(all_covariates .== "intercept_ra3_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra3_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra3_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 3) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra4_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra4_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra4_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra4_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 4) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];

			X[:,findall(all_covariates .== "intercept_ra5_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra5_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra5_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 5) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
					
			X[:,findall(all_covariates .== "intercept_ra6_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 6) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra6_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 6) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra7_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra7_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra7_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra7_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 7) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];	
						
			X[:,findall(all_covariates .== "intercept_ra8_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra8_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra8_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];	
			X[:,findall(all_covariates .== "intercept_ra8_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 8) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra9_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 9) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra9_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 9) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra10_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra10_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra10_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 10) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];		
			
			X[:,findall(all_covariates .== "intercept_ra11_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 11) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra11_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 11) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra12_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 12) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra12_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 12) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra14_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra14_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra14_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 14) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra15_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra15_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra15_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "intercept_ra15_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 15) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra16_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra16_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];							
			X[:,findall(all_covariates .== "intercept_ra16_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "intercept_ra16_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 16) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
				
			X[:,findall(all_covariates .== "intercept_ra17_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra17_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra17_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
			X[:,findall(all_covariates .== "intercept_ra17_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 17) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			
			X[:,findall(all_covariates .== "intercept_ra18_Anthem")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra18_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra18_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];			
			X[:,findall(all_covariates .== "intercept_ra18_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 18) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
				
			#X[:,findall(all_covariates .== "intercept_ra19_Anthem")[1]] =
			#	(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
			#	X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "intercept_ra19_Blue_Shield")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
				X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "intercept_ra19_Kaiser")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
				X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "intercept_ra19_Health_Net")[1]] =
				(1 .- uninsured_plans) .* (households[household_numbers,:rating_area] .== 19) .*
				X[:,findall(all_covariates .== "Health_Net")[1]];
		end
		
		# Residual from reduced form stage
		if control_function
			X[:,findall(all_covariates .== "residual")[1]] = convert(Array,data[!,:residual]);
		end
		
		# Random Parameter Variables
		if random_parameters
			#X[:,findall(all_covariates .== "Anthem_random")[1]] = X[:,findall(all_covariates .== "Anthem")[1]] .* V[:,findall(random_product_covariates .== "Anthem_random")[1],1];
			#X[:,findall(all_covariates .== "Blue_Shield_random")[1]] = X[:,findall(all_covariates .== "Blue_Shield")[1]] .* V[:,findall(random_product_covariates .== "Blue_Shield_random")[1],1];
			#X[:,findall(all_covariates .== "Kaiser_random")[1]] = X[:,findall(all_covariates .== "Kaiser")[1]] .* V[:,findall(random_product_covariates .== "Kaiser_random")[1],1];
			#X[:,findall(all_covariates .== "Health_Net_random")[1]] = X[:,findall(all_covariates .== "Health_Net")[1]] .* V[:,findall(random_product_covariates .== "Health_Net_random")[1],1];
			if !inattention
				X[:,findall(all_covariates .== "intercept_random")[1]] = X[:,findall(all_covariates .== "intercept")[1]];
				#X[:,findall(all_covariates .== "HMO_random")[1]] = X[:,findall(all_covariates .== "HMO")[1]];
			end
			X[:,findall(all_covariates .== "Anthem_random")[1]] = X[:,findall(all_covariates .== "Anthem")[1]];
			X[:,findall(all_covariates .== "Blue_Shield_random")[1]] = X[:,findall(all_covariates .== "Blue_Shield")[1]];
			X[:,findall(all_covariates .== "Kaiser_random")[1]] = X[:,findall(all_covariates .== "Kaiser")[1]];
			X[:,findall(all_covariates .== "Health_Net_random")[1]] = X[:,findall(all_covariates .== "Health_Net")[1]];
			if !add_av_interactions
				X[:,findall(all_covariates .== "AV_random")[1]] = X[:,findall(all_covariates .== "AV")[1]];
			end
			if add_silver
				X[:,findall(all_covariates .== "Silver_random")[1]] = X[:,findall(all_covariates .== "Silver")[1]];
			end
			
			#X[:,findall(all_covariates .== "Premium_random")[1]] = X[:,prem_ind];
			#X[:,findall(all_covariates .== "Premium_random")[1]] = X[:,prem_ind].* V[:,:,1];
		end
		
		return(X)
	end	
	
	#### Term for control function approach	
	if control_function
		eta = Array(randn(length(choices),ns));		
	end	
		
	#### Terms for Random Parameters
	V=0;
	if random_parameters 
		
		# Assign a household reference number
		household_ref_numbers = zeros(maximum(households[!,:household_id]));
		household_ref_numbers[unique(households[!,:household_id])] = collect(1:1:length(unique(households[!,:household_id])));
		households[!,:ref_number] = Int64.(household_ref_numbers[households[!,:household_id]]);
		
		#if add_complex_interactions
		#	V = zeros(length(household_numbers),length(random_variables),ns);
		#	std_normal = MvNormal(zeros(length(random_variables)), Matrix{Float64}(I,length(random_variables),length(random_variables));
		#	Vi = transpose(rand(std_normal, length(unique(households[!,:ref_number]))));
		#	Vi = Vi[households[!,:ref_number],:];
		#	V = Vi[household_numbers,:];
		#else
		
		V = zeros(length(household_numbers),length(random_variables),ns);
		for k in 1:length(random_variables)
			Random.seed!(random_parameter_product_indices[k]);
			Vi = randn(length(unique(households[!,:ref_number])),ns);
			Vi = Vi[households[!,:ref_number],:];
			V[:,k,:] = Vi[household_numbers,:];
		end
		
		if add_av_as_lognormal
			av_index = findall(all_covariates .== "AV")[1];
			av_sd_index = findall(all_covariates .== "AV_random")[1];
			av_random_index = findall(random_variables .== "AV")[1];
			non_av_indices = setdiff(1:1:length(all_covariates),cat(dims=1,av_index,av_sd_index));
		end
	end
	
		
		# V - IJ * K3 * ns
			# This array will be huge, but I'm trying to simplify the computation later
			# The random input should be the same for each household across plans
	
		# Assign a household reference number
		#household_ref_numbers = zeros(maximum(households[!,:household_id]));
		#household_ref_numbers[unique(households[!,:household_id])] = collect(1:1:length(unique(households[!,:household_id])));
		#households[!,:ref_number] = Int64.(household_ref_numbers[households[!,:household_id]]);
		
		#Vi = randn(length(unique(households[!,:ref_number])),length(random_product_covariates),ns);
		#Vi = Vi[households[!,:ref_number],:,:];
		#V = Vi[household_numbers,:,:];
	
			